First-principle studies of the spin-orbit and the Dzyaloshinskii-Moriya 
interactions in the {CU3} single-molecule magnet 



J.F. Nossa, M.F. Islam, and CM. Canali 
School of Computer Science, Physics and Mathematics, Linnaeus University, Kalmar-Sweden 

M.R. Pcdcrson 
Naval Research Laboratory, Washington DC, USA 
(Dated: November 15, 2011) 

Frustrated triangular molecule magnets such as {CU3} are characterized by two degenerate S=l/2 
ground-states with opposite chirality. Recently it has been proposed theoretically [PRL 101, 217201 
(2008)] and verified by ab-initio calculations [PRB 82, 155446 (2010)] that an external electric field 
can efficiently couple these two chiral spin states, even in the absence of spin-orbit interaction (SOI). 
The SOI is nevertheless important, since it introduces a splitting in the ground-state manifold via 
the Dzyaloshinskii-Moriya interaction. In this paper we present a theoretical study of the effect of 
the SOI on the chiral states within spin density functional theory. We employ a recently-introduced 
Hubbard model approach to elucidate the connection between the SOI and the Dzyaloshinskii- 
Moriya interaction. This allows us to express the Dzyaloshinskii-Moriya interaction constant D in 
terms of the microscopic Hubbard model parameters, which we calculate from first-principles. The 
small splitting that we find for the {CU3} chiral state energies (A w 0.02 meV) is consistent with 
experimental results. The Hubbard model approach adopted here also yields a better estimate of 
the isotropic exchange constant than the ones obtained by comparing total energies of different spin 
configurations. The method used here for calculating the DM interaction unmasks its simple funda- 
mental origin which is the off-diagonal spin-orbit interaction between the generally multireference 
vacuum state and single-electron excitations out of those states 



I. INTRODUCTION 

In the last twenty years single-molecule mag- 
nets (SMMs) have been widely studied both for 
their fundamental physical properties^, and for pos- 
sible applications in magnetic storage and quan- 
tum information^^ Unlike traditional bulk mag- 
netic materials, molecular magnetic materials can 
be magnetized in a magnetic field without any in- 
teraction between the individual molecules. This 
magnetization is a property of the molecules them- 
selves. The magnetization occurs because of the 
large ground-state spin and the large easy-axis mag- 
netic anisotropy barrier separating spin-up and the 
spin-down states. In principle it is possible to store 
and manipulate information in one SMM. Further- 
more the two quantum states representing the two 
possible spin orientations can be used to build a 
quantum qubit. Whether used as classical magnetic 
storage units or as quantum coherent elements, the 
crucial requirement in both cases is the ability to 
control and manipulate the magnetic states of the 
SMM in an efficient way. Manipulation by magnetic 
fields is straightforward but, in practice, cannot be 
realized with molecular-size spatial resolution and 
at fast temporal scales. Unlike magnetic fields, elec- 
tric fields are easy to produce, quickly switched and 
can be applied locally at the nano and molecular 
scale. Therefore manipulation of the properties of 



SMMs by external electric fields is an attractive and 
promising alternative^ 

Although electric fields do not directly couple to 
spins, electric manipulation of the spin states is pos- 
sible indirectly via spin-orbit coupling. This requires 
the presence of a strong spin-orbit coupling such 
that the electric field can effectively flip the spin 
states by acting on the the orbital part of the spin- 
orbitals. When SMMs are involved, this is not the 
most efficient mechanism, since the relative strength 
of spin-orbit interaction scales like the volume of the 
molecule. 

Recently, a different mechanism of spin-electric 
coupling in antifcrromagnctic SMMs, characterized 
by lack of inversion symmetry and spin frustration, 
has been proposed A The best example of such a sys- 
tem is a triangular spin s = 1/2 ring with antifcrro- 
magnetic coupling, realized for example in the {CU3} 
SMM. The low energy physics of this system can be 
described by a three-site spin s — 1/2 Heisenberg 
Hamiltonian whose ground-state manifold is com- 
posed of two degenerate (total) spin S = 1/2 dou- 
blets, with wave functions represented by 



\x±,s z = +-) = — (\ ;tt) + c±| Ut> + *fI fU» , 

(i) 



|X±, s, = ~) = -J=(| m) + e±| m) + e T | lit)) , 

(2) 

where the many-body states \a\a20~3) arc products 
of spin-orbital states Oi = (^,|),i = 1,2,3 localized 
on the three magnetic ions of the molecules, and 
e± = exp (±2?ri/3). The four states \x±,S z = ±1/2) 
in Eqs. ([T]), © are labeled by the eigenvalues S z = 
±1/2 of the z-component of the total spin, and by 
the chirality quantum number \± = ±1, that is, the 
eigenvalues of the chiral operator 
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C z = -^=si • s 2 x s 3 . (3) 

An electric field couples to the SMM through 
eE • R, where e is the electron charge and R = 
J2^=i r i- The two spin-orbital states \x±,S z ), cha- 
racterized by opposite chirality and equal spin pro- 
jection, form the basis of a two-dimensional E' irre- 
ducible representation of D 3 h . General group theory 
arguments then guarantee that the matrix elements, 
e(x+i,Sz\ x -\X-i,Sz) = e(x-i,S z \X + \x+i,S z ) = 
2id 0, where X± = ±X + iY are the in-plane 
components of R, which also transform as the two- 
dimensional irreducible representation E'. Here d 
is a real number that is refereed to as spin electric- 
dipole coupling. It follows that, due to these non- 
zero matrix elements, an electric field can cause tran- 
sitions between two ground state wavcfunctions of 
opposite chirality but with same S z . 

The observation of such electric-field induced 
transitions from one chiral state to another requires 
that the degeneracies between these states be lifted. 
The anisotropic Dzyaloshinskii-Moriya (DM) inter- 
action plays a crucial role in that, it provides one 
possible mechanism that lifts the degeneracies be- 
tween states of different chirality without mixing 
them, as shown in Fig. [TJ More in general, the 
presence of DM interaction provides a mechanism 
to control the size of quantum entanglement in mag- 
netic trimers as a function of the temperature and 
external magnetic field^. Experimentally the DM- 
induccd splitting in {CU3} is estimated to be small 
(approximately 0.5 K— ). 

Recently^ we have investigated the details of the 
electronic properties of the {CU3} SMM, which 
is one of the most promising triangular spin 1/2 
molecules where the spin-electric effect can be rea- 
lized. In particular, we introduced a scheme to 
evaluate the strength of spin electric-dipole cou- 
pling d using ab-initio methods. However, the value 
of the anisotropic DM-exchange constant interac- 
tion, which is responsible of the GS zero-field split- 
ting, has not yet been calculated. The purpose of 
this work is to calculate this splitting by ab-initio 
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Magnetic field 



FIG. 1: (Color online) Schematic diagram of electric- 
field-induced transitions between states of different chi- 
rality belonging to the spin S = 1/2 ground-state mani- 
fold of a triangular antiferromagnet. A represents zero- 
field splitting of the chiral states due to Dzyaloshinskii- 
Moriya interaction. 



methods. In order to achieve this goal, we ana- 
lyze the microscopic origin of the DM interaction 
via a Hubbard model approach in the presence of 
spin-orbit integration, which is the correct mini- 
mal model to describe both spin and charge fluc- 
tuations of these strongly correlated electron sys- 
tems. At half- filling and in the large Hubbard U 
limit, spin-dependent virtual hopping processes, in- 
duced by the spin-orbit interaction, give rise to 
an anisotropic exchange interaction.— There is a 
close analogy with the isotropic Hciscnbcrg exchange 
interaction obtained in second-order perturbation 
theory in the spin-independent hopping perturba- 
tion. Beside elucidating the physical mechanism 
leading to the anisotropic DM exchange interaction, 
this approach provides a very convenient prescrip- 
tion on how to extract the DM exchange constant 
from first-principle calculations, which we have ca- 
rried out for {CU3}. 



This paper is organized as follows. In Sec. Ill Al we 
discuss the general properties of the DM interaction. 
The Hubbard model approach for calculating DM 
vector, adopted in this work, is discussed in Sec. Ill Bl 
In Sec. IIIII wc discuss details of extracting Hubbard 
model parameters from our ab-initio calculations. In 
Sec. MI PI we discuss other methods that are usually 
employed for calculating the DM vector. Finally in 
Sec. IIVI we present a summary of our work. 
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II. THE DZYALOSHINSKII-MORIYA 

INTERACTION IN FRUSTRATED 
ANTIFERROMAGNETIC SPIN RINGS 

A. General properties of the DM interaction 

The Dzyaloshinskii-Moriya (DM) interaction is an 
anisotropic exchange interaction resulting from the 
interplay of the Coulomb interaction and the spin- 
orbit coupling in systems of low crystal symmetry. 
The DM interaction is an important effect for many 
magnetic systems and plays a crucial role in deter- 
mining the zero-field splitting of energy levels. An 
anisotropic exchange interaction of the form 

D 12 • Si x S 2 , (4) 

which is linear in the spin-orbit interaction, was 
first put forward by Dzyaloshinskii on the basis of 
symmetry considerations^ Later Moriy a 10 ' 11 pro- 
vided a mechanism for this interaction by exten- 
ding Anderson's theory of superexchangoi^ to in- 
clude the effect of spin-orbit coupling. Let us con- 
sider for simplicity two "magnetic ions", R and R', 
each occupied by a single electron in the ground 
state. Second-order perturbation theory in the ho- 
pping Hamiltonian H t coupling the two sites gives 
rise to an isotropic antiferromagnetic interaction 
with exchange constant J = 2t RR ,/U, where t RR < 
is a spin-independent hopping integral and U is 
the energy required to transfer an electron from 
R to R'. When spin-orbit interaction -ffsoi is in- 
cluded, similar second-order processes can generate 
an anisotropic exchange interaction in the form of 
Eq. (g]), with D ~ t RR ,b RR ,/U where 6rr- is a 
SOI-induccd or SOI-dcpendcnt hopping integral. To 
lowest-order, 6rr/ is just the matrix elements of the 
Hsoi between two orbitals localized at R and R'. 
This is the dominant contribution to D. In case 
that at each site more than one orbital |R, /i) , /i = 
1,2,- •• plays a role, higher-order terms such as 
b RR > = t RR , (R, n\H SO i\R' , v') / AE^ are possible, 
making the corresponding D effectively a third-order 
coupling in the perturbations H t and -ffsoi- It turns 
out that D ~ (Ag/g), where g is the free-electron 
gyromagnetic ratio and Ag the deviation from g in- 
duced by SOIJJ. 

As shown by Moriya, other terms linear in the 
SOI contribute to the anisotropic exchange of the 
form of Eq. The second most important con- 
tribution is also a second-order term resulting from 
SOI and direct inter-atomic exchange interaction 
J CX (R, R'). In antiferromagnetic crystals this term 
is J CX (R, R')/J times smaller than the second-order 
contribution proportional to t RR >b RR i . Finally, 
third-order contributions to D include the hopping 



terms twice and the intra-atomic exchange constant 
Jo- They are Jq/U smaller than second-order terms. 

The DM exchange vector D vanishes when the 
symmetry of the crystal is high. This is the case, for 
example, when the point located halfway between 
the two magnetic ions in a unit cell is a center of 
inversion. In low-dimensional crystals where D^O, 
the anisotropic exchange is typically the most im- 
portant anisotropic contribution between spins. The 
DM interaction favors non-collinear spin configura- 
tions, with typical canted spins. As such, it de- 
termines the spin arrangements and it is responsi- 
ble for the weak ferromagnetism observed in some 
predominantly antiferromagnetic crystals such as a- 
ES2O2. The tendency toward canted spin configura- 
tions is most-easily seen by minimizing the energy in 
Eq. (|4]) for two classical spins, when the DM vector 
D is, for example, perpendicular to the line joining 
the two ions. It can be shown that the minimum 
energy corresponds to a spin configuration where 
both spin are perpendicular to each other and to 
the direction of D. Similar conclusions can be ob- 
tained by analyzing the same system quantum me- 
chanically. The DM interactions is also responsi- 
ble for proposed non-collinear spin configurations in 
magnetic clusters engineered by STM techniques on 
insulating surfaces ! 13 ' 14 

B. The DM interaction for antiferromagnetic 
spin rings within a Hubbard model approach 

In this section we specialize the previous discus- 
sion to the case of an antiferromagnetic spin triangle, 
and show how the DM interaction can be derived mi- 
croscopically from a Hubbard model at half filling, 
in the presence of spin-orbit interaction. 

As mentioned in the introduction, the low-energy 
magnetic properties of {CU3} are well-described by 
an isotropic antiferromagnetic Heisenberg model 

H R = JijSi ■ sj , Jij > , (5) 

where Sj are spin vector operators of magnitude 
Si = 1/2, predominately localized at the three Cu 
sites. If the small distortion from a perfect equilat- 
eral arrangement of the three Cu atoms is neglected, 
the three exchange constants are the same, Jij = J. 
DFT calculations! find J w 3.7 meV. The GS man- 
ifold comprises two spin 5 = 1/2 doublets, which 
can be represented by the two chiral states given 
in Eqs. (flj, ([2]), or any two orthogonal linear com- 
bination of these. The spin 5 = 3/2 excited-state 
multiplet is separated by the GS by an energy of 
order J. 
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It is well-known that the AFM Heisenberg model 
represents an effective low-energy spin model that 
can be derived from an underlying Hubbard model 
at half-filling in the large t/U limit. The choice 
of the best minimal model capturing the essential 
microscopic features of the electronic system is of- 
ten a complex task, particularly when the exchange 
interaction between the magnetic ions is mediated 
via several paths involving non-magnetic ions, as for 
the case of {CU3}. We will neglect these complica- 
tions and assume that an effective one-band Hub- 
bard model suffices for this purpose. We will see 
that our first-principles calculations corroborate this 
choice, showing that one localized orbital at each 
magnetic ion indeed is enough to describe the low 
energy physics of the system. We will comment 
later on the possibility of considering a more com- 
plex Hubbard model to describe the non-magnetic 
bridges between Cu atoms, as well as the need of 
including more than one orbital at the Cu sites. 

The second quantized one-band Hubbard Hamil- 
tonian reads 

H v = -t E E {4*S'a + h - c -} + \ U J2 n ^ n ^ ' 

i.j a i 

(6) 

where c\ a (c\ a ) creates (destroys) an electron with 
spin a at site i and rii a = c\ a Ci a is the particle 
number operator. More precisely the index i la- 
bels a Wannier function localized at site i. The first 
term represents the kinetic energy, characterized by 
a spin-independent hopping parameter t, which is 
the same for all pairs of site due to the C3 symme- 
try of the CU3 molecule magnet. The second-term is 
an on-site repulsion energy of strength U, which has 
an effect only when two electrons of opposite spins 
reside on the same site. It is the on-site repulsion 
energy. 

The spin-orbit interaction in the Hubbard model 
is described by adding the following spin-dependent 
hopping term2»iSr— 

^01 = E E {4 ir§- ■ <><# ) c iP + h - c ] > ( 7 ) 

where a is the vector of the three Pauli matrices. 
Here the vector Pjj is proportional to the matrix 
element of V7 x p between the orbital parts of the 
Wannier functions at sites i and j; V is the one- 
electron potential and p is the momentum operator. 
Clearly the spin-orbit term has the form of a spin- 
dependent hopping, which is added to the usual spin- 
independent hopping proportional to t. This form of 
the spin-orbit interaction is a special case of Moriya's 
hopping terms^ in the limit that all but one orbital 
energy is taken to infinity.— 



In contrast to the spin-independent hopping term, 
the spin-depending hopping parameters are related 
by both the full symmetry of the molecule and the 
local symmetry of localized orbitals£ Now, because 
of the o~ v symmetry, = Pe z . The final expres- 
sion of the Hubbard model, including the spin-orbit 
interaction is 

Hu+soi = E { c «*( ~ 1 + * A soia) c l+ ia + h.c.j 

+2 U ^2 n it n a > ( 8 ) 

i 

where Asoi = P/2 = Py/2 • e z is the spin-orbit 
parameter. 

We want to treat the two hopping terms pertur- 
batively on the same footing, by doing an expansion 
around the atomic limit t/U , \§oi/U — > 0. In many 
molecular magnets t 3> Asoi- This turns out to be 
the case also for {CU3}. In other molecules the two 
hopping parameters are of the same order of magni- 
tude. 

We are interested in the half-filling regime. We 
know that second-order perturbation theory in t 
results in an antiferromagnetic isotropic exchange 
term that splits the spin degeneracy of the low- 
energy sector of the Hubbard model, defined by 
the singly-occupied states. This action can be rep- 
resented with an effective spin Hamiltonian, the 
isotropic Heisenberg model, with exchange constant 
J = 4:t 2 /\U\*£ Similarly Loss et al. showed that 
another second-order term proportional to tXsoi/U 
generates an anisotropic exchange term that can be 
identified with the DM interaction^ They write ap- 
proximate adapted many-body states to first-order 
in the perturbation \t\, Asoi "C U, corresponding to 
singly-occupied states. In particular there are two 
independent doublets, 

Whi) = ^=(\ m)+e±\ m) + e T \ m» , (9) 

and 

ivi?) = ^=(i m> + e±i m) + <m ut)) , (io) 

with e± = exp (±27ri/3). These are states with 
5 = 1/2 and S z = ±1/2. These states are formally 
identical to the chiral states given in Eqs. ((J) and 
^ . Now, each of the terms appearing in these equa- 
tions is a single Slater determinant obtained by three 
creation operators acting on the vacuum, e.g., 

I tU) = clt4f4j > ■ (11) 
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The states \i(>e?) and are eigenstates of Hub- 

bard Hamiltonian when t = Agoi = 0. The tunnel- 
ing and SOI mix the singly-occupied and double- 
occupied states. The first-order correction is ob- 
tained by mixing in doubly-occupied states 



| (e--l)(*±«A SO i) ^ 



V2U 



2a 
E' 1 



3e + (t ± aXsoi) 
V2U 



(12) 



and allows us, in the low-energy regime, to make the 
identification 



5A SO it 
U 



(18) 



This Hubbard model analysis suggests an avenue 
to extract the DM parameters from an ab-initio cal- 
culation. Only three parameters are needed, namely 
the spin-orbit interaction Xsoi > the hopping param- 
eter t and the on-site repulsion energy U . 



where 



2a 
B'l 



V6 



3 



1.2 



Hm+m)) , as) 



and 



V i 



(14) 



with = c} t c},c]JO) (i = 1,2,3 and j £ i) rep- 
resenting the double-occupied sites. 

The next step is to take the expectation value of 
the spin-orbit part of Eq. (JSJ in these approximated 
states. The result is£ 



1 a i 



± s gn(a) ■ (15) 



Note that off-diagonal matrix elements of -ffsoi 
vanish; in other words, SOI splits but does not mix 
the chiral states. 

In the small t/U , Xsoi I U limit, we can resort to a 
spin-only description of the low-energy physics of the 
system. The ground-state manifold (corresponding 
to the states Eq. (|12[) ) is given by the two chiral spin 
states Eqs. ©, 0. 

The anisotropic DM spin exchange Hamiltonian 
in I?3h symmetry is given by^ 



DM 



iD z yr^ 
i 



i+l i 

S_ — S_ 



(16) 



Now, in the low energy regime corresponding to 
a £*3h symmetric molecule magnet, the spin-orbit 
interaction can be reduced to the effective form 



H, 



SOI 



= A 



SOI 



C Z S Z 



(17) 



where A$oi is the effective SOI coupling constant. 

The DM interaction expressed in this form clearly 
shows that it splits but does not mix the two chiral 
states.— The splitting is exactly proportional to D z 



C. Semiclassical analysis of the DM 
interaction in frustrated spin systems 

The quantum mechanical frustration present in an 
antiferromagnetic spin triangle and the DM interac- 
tion both tend to favor non-collinear spin configu- 
rations. It is instructive to study their interplay in 
a semiclassical approach, where non-collincarity is a 
more intuitive concept. 

The classical Heiscnberg model with energy func- 
tional given by Eq. ((5|) has two degenerate "ground 
states" , given by the two non-collinear spin configu- 
rations shown in Fig. [2] Classically these two states 
are the best way to by-pass the frustration present 
for any collinear spin configuration in a triangular 
antifcrromagnct. Quantum mechanically the two 
non-collinear spin configurations can be represented 
by the states 



nc±; 



(ai|t)i+/3iU)i)®(a2|t)2+^U> 2 ) 



)(a 3 |t>3 + /%U>3) 



(19) 



where a = cos(#/2) and (3 = exp{i0} sin(0/2). Here 
9 is the elevation angle and is the azimuth angle. 
The three spinors (a»| f)i + Pi\ , i = 1,2,3 are 
three spin- 1/2 coherent states defined by the three 
non-collinear directions obtained by rotating consec- 
utively by the angle ±240° (see Fig. [2j). Anticlock- 
wise rotations (by —240°) define a left-handed heli- 
cal state (Fig. 2(a) |; clockwise rotatio ns (by +240°) 



2(b) I 



define a right-hand helical state (Fig. 

In contrast to the true GS given in Eqs. (QJ) , (0) the 
non-collinear states defined in Eq. (fT"9)) are neither 
eigenstates of the quantum Hamiltonian Eq. ([5]) nor 
of S 2 and S z . The expectation value of the Hamil- 
tonian Hh at these states is defined by 



(VVc±I-ffH|V>: 



nc±; 



3 J/4 . 



(20) 



The fact that the energy of the collinear states is 
higher than the energy of the chiral states by 3 J/8 is 
not surprising, since the noncollincar states defined 




III. AB-INITIO CALCULATION OF THE 
DM VECTOR 




(a)Left-handcd 



(b)Right-handed 



FIG. 2: Non-collinear helical system of three Cu-atom 
spins. Anticlockwise rotations (by —240°) define a left- 
handed helical state (a); clockwise rotations (by +240°) 
define a right-hand helical state (b). 



in Eq. (QIJJ) arc a mixture of S = 1/2 and S = 3/2 
components. 

When rewritten in term of the electronic states for 
the corresponding Hubbard model at half-filling in 
the small t/U limit, the non-collinear spin-coherent 
states defined in Eq. (|19|) can be considered to be the 
"best" energy states given by a single Slater determi- 
nant (Note that the chiral states cannot be written 
as a single Slater determinant). 

It is now interesting to examine the effect of 
the DM interaction on these states. A straightfor- 
ward calculations shows that for the DM interaction 
of Eq. (|16[) . where only the z-componcnt of D is 
nonzero 



(V'nc±|#DM|l/'nc±! 



= ± 3 -^D z 
4 2 



(21) 



Therefore, as for the GS manifold of the ex- 
act eigenstates, the DM interaction splits but does 
not couple the two noncollinear states. The DM 
parameter D z is, by Eq. (f2"Tj) . related to the 
DM interaction-induced energy-gap between the two 
noncollinear states 



AE nc = (ip nc + |#dm|VVc - 



nc — / 



(22) 



This result suggests a way of extracting the DM 
vector parameter D similar in spirit to the method 
used to calculate the isotropic exchange parameter 
J by comparing the energy difference of states with 
ferromagnetic and antiferromagnetic spin configura- 
tions respectively. In the next section we will see 
that this procedure can also be carried out by first- 
principle methods. 



All the calculations in this work are carried out by 
using ab initio package NRLMO L 19 ' 20 , which uses a 
Gaussian basis set to solve the Kohn-Sham equa- 
tions within PBE-GGA approximation.— For more 
computational details and the electronic properties 
of {CU3} we refer the reader to our previous work^ 



A. Calculation of the hopping term t 



As discussed in the section III B[ the Hubbard 
model approach is based on allowing the localized 
electrons to hop to its nearest neighbor sites and in 
the present case of the {CU3} molecule, these local- 
ized electrons are d electrons. Therefore, for cal- 
culating hopping parameter t, the relevant states 
are those d electron states that lie close to the 
Fermi level. Let \K, a) be the three relevant Kohn- 
Sham eigenstates calculated from NRLMOL. We can 
write them as a linear combination of the local- 
ized atomic orbitals, centered at the three Cu sites, 
{\<j> a ) , \4> b ) , \(j) c )} (g) \xa), with a =t,i for spin up 
and down, respectively: 



\K,a) 



(23) 



where C % Ka is the weight of the localized |<fo) \x a ) 
wavefunction. 

For the |ttt) s P m configuration, in the absence 
of spin-orbit interaction, the relevant three levels 
around the Fermi level are doubly and singly de- 
generate. These levels are sketched in Fig. [3] 



Fermi 
energy 



\A,i) 



FIG. 3: Schematic diagram of the Kohn-Sham energy 
levels around the Fermi level 

We obtain the level structure by diagonalizing the 
three-site Hamiltonian in the absence of the SOI: 

H =e J2\<t>i)(<f>i\-tJ2 l&X&l , (24) 
i i=£j 

where Eq is the onsite energy, t is the hopping term 
and i,j = a, 6, c represent the copper sites. We get 
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the eigenvalues So + t and Eq — 2t for the two-fold and 
one-fold degenerate states, respectively. The Kohn- 
Sham eigenvectors can be defined as a linear combi- 
nation of the localized wavcf unctions, 

\ElA) = (|0a> - \<M) It) , 

1^2, t) = -^(l0«> + l^>-2|^ c ))|t> , (25) 

I At) = (|*«> + \h) + \<f>c)) It) • 

Now the localized states can be written in term of 
the Kohn-Sham functions 



\<t>a) It) 

V J V ^ V o 

l«l t) = ^t)_^) + ^t>, (26) 

l« It) 



lAt) 


+ 1 Y } + 


\E 2 ,t) 






V6 


lAt) 




\E*A) 


73 






lAt) 






V3 


V6 





Our calculations showed that these states are pri- 
marily localized on the Cu atoms and have d charac- 
ter. We have obtained the Kohn-Sham eigenenergics 
for the one-fold and two-fold degenerate states 

(£ ljT |ffo|£l,t) = \((<t>a\-(<Pb\)H (\cf>a)-\^)) 

= £a + t , 

(A,t\H \A,t) = i((^| + (^| + (^|)Ho 

= e -2t. (27) 

From Eqs. (f2"7} we can finally evaluate the value 
of the parameter t as: 

* = i((£i,tl#o|£i,t>-(At|ffo|At)) 
= 50.84meV. (28) 



B. Calculation of the spin-orbit interaction 
parameter Xsoi 

Standard spin-orbit interaction representation for 
spherical systems is given by 

CUr,L,S) = — S-L--A1, (29) 

where r is the position, L is the angular momentum, 
S is the spin moment, c is the speed of light, and $ is 
a spherically symmetric potential. The above equa- 
tion is exact for spherical systems. For a multicen- 
ter system a superposition of such terms needs to be 



considered. However, this approximation could miss 
non-spherical correlations important for anisotropic 
energies. Instead of using Eq. (|2^|) . a generaliza- 
tion of the spin-orbit interaction for non-spherical 
or multicenter systems is given by 

[/ so (r,p,S) = -is.pxV$(r), (30) 

where p is the momentum operator and a external 
electric field is given by E = — V<&. 

Pedcrson et. al (see Ref. l22t ) have shown an exact 
simplified method for incorporating spin-orbit cou- 
pling into density-functional calculations. In order 
to get the basis-set for the spin-orbit coupling the 
single-electron wave function can be expressed as 



^ s (r) = ^C-/ J (r) Xa , 



(31) 



where fj (r) is a spatial basis function, \a is either a 
majority or minority spin spinor, and Cj S a are deter- 
mined by effectively diagonalizing the Hamiltonian 
matrix. In order to calculate the effect of the SOI 
(Eq. (|3"0"|) ) it is necessary to calculate matrix ele- 
ments of the form 

U ja ,ka' = (fjXa\U(T,p,S)\fkXa') 

= 'E,Uf S \V g \f k )(x a \S x \ Xa >) ,(32) 



where 

(33) 

The matrix elements for V y and V z are obtained 



by cyclical permutations of x,y and z in Eq. (|33p . 
This methodology for the SOI matrix gives several 
advantages, namely, it does not require the deter- 
mination of the electric field; it is specially ideal for 
basis functions constructed from Gaussian-type or- 
bitals, Slater-type functions, and plane waves. 

We are interested in the matrix elements in the 
localized basis-set, Eq. (|26|) : 

(^l(Xtl^cl^)lxt) = -^5<&|pxV*(r)|0 fc ) 

• (xtl s Ixt) 



2i 



{<t>i\V z \<j> k 



= -^Pk = -i^soi 



(34) 



We can write these matrix elements in the Kohn- 
Sham basis set 

(^l(xtl^ol^)lxt) = E(^t)*^'t 

KK' 

x (A',t|(7so/|^',t)(35) 
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We have obtained the matrix elements for 
the spin-orbit interaction in the Kohn-Sham ba- 
sis, {{Eh) ,\En) ,\A)} <g> \x a ) (Eq. dSll)), and used 
Eqs. (Jini) to obtain the matrix elements: 



0.85 0.85^ 
0.85 
0.85 



(36) 



From Eq. (O we have A 50/ = p z ik /2 = 0.43 meV. 



C. Calculation of the Hubbard U and 
evaluation of D z and J 

The most common approach for calculating U in- 
volves calculation of energy, E, of the molecule with 
N, N + 1 and N — 1 electron and extracting U from 
the equation below, 

U = E(N + 1)+E(N -1)-2E(N) 

= [E(N + 1) - E(N)} - [E(N) - E(N - 1)] 
= A — I . (37) 

In the above equation A is (minus) the electron 
affinity 2 ^ and I is the ionization energy. For systems 
that are not closed shell, such as those considered 
here, the U value is essentially the second derivative 
of energy with respect to charge and it is possible to 
determine U by calculating the energy as a function 
of charge. 

For the single-band Hubbard-model correspond- 
ing to the {CU3} molecule, we are interested in ob- 
taining energies for the charge-transfer excitations 
involving the transfer of a localized rf-electron on 
one copper site to a localized d-electron on another 
site. Specifically we wish to know the energy of 
I X) = |t<4otc) relative to |ta-Utc)- There arc a to- 
tal of twelve charge-transfer excitations that can be 
made with one-site doubly occupied and one elec- 
tron on one of the other sites. For the half-filled 
case of interest here, the energy difference depends 
upon the electron affinity of the state on site a, 
the ionization energy of the state on site b and the 
residual long-range coulomb interaction between the 
negatively charged electron added to site a and the 
positively charged hole that is left behind on site 
b. Since site b and site a are equivalent, it fol- 
lows that we simply need to calculate U for any 
one of the copper sites in the half filled case. A 
very rough estimate of the charge transfer energy 
may be determined by calculating the PBE-GGA 
energy of the Cu atom with an electron configura- 
tion of ls 2 2s 2 3s 2 4s 2 p 6 3p 6 3d™ with n=8,9,10. Using 
n=9 as the reference state, one finds a bare U value 



of 13.76 eV which, after accounting for the particle- 
hole interaction (27.2116/i?cu-c« = 2.95eV, where 
Rcu-Cu — 4.87 Bohr is the distance between mag- 
netic centers) is shifted to 10.8 cV. 

In the {CU3} molecule, we have chosen to cal- 
culate U quasi-analytically by gradually adding (or 
subtracting) a small fraction of electronic charge Sq 
to one of the half-filled Cu d-states. The energy of 
the system as a function of Sq is shown in Fig. [4j 
where we can see that it can well be reproduced by 
a quadratic fitting curve. The figure shows that, 
upon adding a fractional charge to a localized or- 
bital, the total energy initially decreases, since the 
orbital energy is negative. Eventually, however, the 
competing Coulomb repulsion takes over and the net 
change in total energy for adding one electron to a 
localized orbital is positive. In contrast, with one 
extra electron delocalized throughout the molecule 
the total energy is usually smaller than the energy 
of the neutral molecule. 
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FIG. 4: (Color online) Dependence of the total energy 
on added fractional charge Sq. The (blue) circle repre- 
sent the results of NRLMOL calculations and the dashed 
(red) line represents a quadratic fit. 



The difference in the energy of the system before 
and after adding a fraction of electronic charge Sq is 
given by AE = U eff = USq 2 - e 2 Sq 2 /i?c u -Cu, where 
U = d 2 E(q)/dq 2 . Wc have calculated the effective 
parameter U e ff by setting Sq = 1: 



U ef f = Sq 



d 2 E{q) 



dq 2 -Rcu-cu , 
= 9.06 eV , (38) 

where E{q) = E + {U/2)(q - q Q ) 2 with E being a 
constant. 
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1 . Evaluation of D z and J 

Having calculated the parameters t, Aso, and 
U e ff, we are now able to use Eq. (TTg)) and evalu- 
ate the Dzyaloshinskii-Moriya parameter D z . We 
obtain 



example, the states associated with the system de- 
picted in Fig 2(a) would be represented according to: 



D z 



' C/cff 



O.OlmcV 



(39) 



This value of D z yields a small splitting of the chi- 
ral state, A rj 0.02 meV rj 0.3 K. Experimental esti- 
mates of the DM parameter find a splitting 3-4 times 
larger than this value. Considering the smallness of 
this energy and the uncertainty in the experimen- 
tal measurements, the two estimates are consistent 
with each other. On the other hand, it is also possi- 
ble that part of the discrepancy between theory and 
experiment is due to the fact that other mechanisms, 
different from the DM interaction, contribute to the 
splitting. In particular Ref<^ pointed out that small 
deformations of the triangular molecule can lift the 
chiral degeneracy and this contribution to the split- 
ting could be even more important than the DM 
interaction. If this is indeed the case, our results 
would imply that our method of computing the DM 
parameter is actually rather accurate. 

From a computational point of view, it is interest- 
ing at this point to evaluate the isotropic exchange 
constant J from the Hubbard model pcrturbativc 
approach, which gives 



J = At 2 /U lmeV. 



(40) 



This estimate of J is considerably closer to the 
experimental value of 0.5 meV than the value of 3.7 
meV obtained by computing the energy difference 
between states with ferromagnetic and antifcrromag- 
netic spin configurations.— 



D. Comparison with other methods 

In a recent work Takcda et. al£^ have used a 
non-collinear approach to estimate the DM interac- 
tion. Instead of the use of simple product functions, 
this work capitalizes on the use of generalized 
orbitals which are composed of a linear combination 
of both spinors with different and variable spatial 
functions. By using such a representation it is 
possible to develop single-determinants which are 
composed of a linear combination of the chiral spin 
1/2 states and the non-chiral spin 3/2 states. For 



\xix h ± x c ± ) 
1 



lttt>±*i4tt>T(-i) 1/6 mt) 



2V2 

T (-i) 5/6 |ttl) T i 1444) - 1444) 
+(-i) 1 / 3 im) - (-i) 2/3 1441)1 (4i) 



where X+(0, 4>) = cos(0/2) |f) + exp{i0} sin(0/2) \i) 
and X_(0,0) = sin(0/2) |t) - exp{i</>} cos(0/2) ||), 
with = tt/2 and (j) = tt/2, 7tt/2, -tt/2. They fur- 
ther claim that AE nc = i^/AD z (see Eq. |52J) can 
be estimated by a perturbational treatment of the 
SOI, as follows 

A£" nc = (VVc + |-HsOl|V'nc+) - (V'nc-l-H'sOllV'nc-) , 

(42) 

where -ffgoi is the one-electron spin-orbit interac- 
tion. These expectation values can be calculated by 
DFT. 

It is clear from the expression Eq. ([4"Tj) that the 
expectation value of the spin-orbit interaction for 
this and other states would be linear so, without 
other considerations, one can not extract an inter- 
action that depends upon the excitations of interest 
to the Hubbard Hamiltonian. However, in analogy 
to the expansion of the many-electron wavefunction 
for molecular hydrogen in regions intermediate be- 
tween the bonding and separated-atom limit, a self- 
consistent optimization of such a starting determi- 
nant allows the spin-orbitals to be intermediate be- 
tween the doubly occupied and single occupied rep- 
resentations. While the resulting noncollincar wave- 
function is still a single Slater determinant in char- 
acter, expansion of the noncollinear state in terms 
of the Hubbard states would show a wavefunction 
comprised primarily of the 8x8 half-filled determi- 
nants but would also contain small contributions of 
the ionic contributions which are shifted upward by 
U e ff. It is the small admixture of these states that 
allow Takeda et. al. to extract both the exchange 
parameters and the DM interaction through the use 
of noncollinear representations. This approach could 
have advantages from an operational viewpoint since 
it effectively addresses the potential role of other 
excited states that are routinely excluded from the 
Hubbard Hamiltonian. However, the precise inter- 
actions which ultimately mediate the appearance of 
the DM interaction require additional analysis which 
is every bit as arduous as that presented here. 

An alternative method to calculate the DM vec- 
tor, based on Andersen's "local force theorem"—, 
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was developed by Solovyev et al££- More recently 
this method was utilized in conjunction with DFT 
to study the DM interaction between magnetic 
atoms inserted in different crystalline systems and 
surfaces ! 14 ' 27 Essentially this method expresses the 
DM vector in terms of the Green's functions of the 
system, modified by the spin-orbit interaction. Al- 
though computationally sophisticated, the Green's 
function method is physically less transparent than 
the one adopted here, particularly for a finite system 
such a triangular SMM, where the crucial ingredi- 
ents leading to the anisotropic DM exchange can be 
reduced to a few parameters that have a direct phys- 
ical interpretation within the Hubbard model. 

IV. CONCLUSIONS 

We carried out a first-principles investigation of 
the zero-field splitting of the chiral ground states 
of a {CU3} single-molecule magnet (SMM), caused 
by the Dzyaloshinskii-Moriya interaction. Our ap- 
proach relies on the perturbative analysis of a Hub- 
bard model, which includes spin-orbit interaction. 
In the large U limit, appropriate for {CU3}, it is 
possible to express the Dzyaloshinskii-Moriya con- 
stant in terms of the parameters that define the 
Hubbard model, such as the effective hopping inte- 
gral between magnetic sites t, the on-site repulsion 
energy U, and the strength of the spin-orbit Asoi- 
We then carried out an approximate method to ex- 
tract the values of these parameters from our spin 
density functional theory calculations of the SMM. 
The value of the Dzyaloshinskii-Moriya constant D 
that we found is of the order of 0.01 mcV, which is 
a factor of 5 smaller than the value measured ex- 
perimentally. Given the uncertainty of the experi- 
mental result and the fact that other effects might 
contribute to the zero- field spin splitting of the chiral 



states, our estimate should be considered consistent 
with experiment. 

The method of computing the DM parameter by 
effectively extending Anderson's theory of superex- 
change to include spin-orbit interaction is very close 
to Moriya's original formulation of anisotropic ex- 
change. It is interesting to note that if we use 
this approach to calculate the isotropic superex- 
change constant J of the Heisenberg model describ- 
ing {CU3}, we obtain a value that is closer to exper- 
imental result than the estimates based on total en- 
ergy calculations of ferromagnetic vs antifcrromag- 
netic spin configurations. This seems to suggest that 
this approach is not only physically very intuitive, 
but it might also bear promise of good numerical 
accuracy. 

While the methods discussed here provide physi- 
cal insight into the nature of the DM interaction, we 
note that for future calculations it would be desir- 
able to consider excitations that are not normally 
included in the single-band Hubbard model. For 
such an approach it would be necessary to include 
methodologies that allow for the calculation of all 
excitations in such systems. 
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